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A METHO D OF RADAR PATTERN RECOGIMITIQN BY SORTING 
SIGNALS INTO DATA CLUSTERS 

TECHNICAL FTFI D 

The present invention relates, in general, to radar identification and, 
more specifically, to sorting-signals received from multiple emitters into data clusters 
for pattern recognition. 

BACKGROU ND OF THE INVENTION 

Radars emit a variety of signals that may characterize and identify 
them. Each radar may emit a specific pulse amplitude and a specific fixed radio 
frequency (RF) or a variable RF ranging over a fixed bandwidth. Each may emit a 
fixed pulse repetition frequency (PRE) or pulse repetition interval (PRI) and may be 
of a certain pulse width (PW). 

An aircraft flying into a region with an onboard wideband receiver may 
detect a variety of signals emitted from multiple radars located in that region. 
Unless these signals are sorted and separated from each other, it is not possible for 
the aircraft to determine the types of classes of radars it is about to encounter. It 
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does not know whether the radars are hostile and does not know whether the radars 
present a high or low threat to the incoming aircraft. 

A need, therefore, exists for an aircraft to be able to sort and identify 
the variety of radars that are emitting energy towards the aircraft. The present 
invention addresses this need. 

SUMI^IARY OF THE INVEIMTTniM 

To meet this and other needs, and in view of its purposes, the present 
invention provides a method of classifying radar emitters including the steps of: (a) 
receiving a plurality of signals from the radar emitters; (b) generating data 
components for each signal received from the radar emitters; (c) forming multi- 
dimensional samples using the generated data components; and (d) sorting the 
multi-dimensional samples into a plurality of data clusters, based on their respective 
proximity to the data clusters, each data cluster representing a classification of a 
radar emitter. Step (b) includies generating pulse data descriptors (PDWs) during a 
predetermined interval of time, and generating at least radio frequency (RF) data 
and pulse width (PW) data for the PDWs. 

In another embodiment, the invention provides a system for classifying 
radar emitters including a receiver for receiving a plurality of signals from the radar 
emitters, and a processor coupled to the receiver, for (a) generating data 
components for each signal received from the radar emitters, (b) forming multi- 
dimensional samples from the generated data components; and (c) sorting the multi- 
dimensional samples into a plurality of data clusters, based on their respective 
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proximity to the data clusters, each data cluster representing a classification of a 
radar emitter. 

The processor may also assign a multi-dimensional sample to a data 
cluster, based on a Euclidean distance between the multi-dimensional sample and a 
center of the data cluster. The center of the data cluster may be formed as a mean 
vector of a set of multi-dimensional samples assigned to the data cluster. 

In yet another embodiment, the invention provides a machine readable 
storage medium containing a set of instructions for a computer. The set of 
instructions implements the following steps: (a) processing a plurality of signals 
received from a receiver; (b) generating data components for each signal received 
from the receiver; (c) forming multi-dimensional samples using the generated data 
components; and (d) sorting the multi-dimensional samples into a plurality of data 
clusters, based oh their respective proximity to the data clusters, each data cluster 
representing a classification of a radar emitter. Step (d) may include sorting the 
multi-dimensional samples using an ISODATA (iterative self-organizing data analysis 
technique) computer algorithm. 

It is understood that the foregoing general description and the 
following detailed description are exemplary, but are not restrictive, of the invention. 
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BRIEF DESCRIPTION OF THE DRAWING 

The invention is best understood from the following detailed 
description when read in connection with the accompanying drawing. Included in the 
drawing are the following figures: 

FIG. 1 Is a flow diagram showing steps performed by a conventional 
ISODATA computer algorithm; 

FIG. 2 is a system block diagram illustrating a radar threat sorting, 
clustering and classification system, in accordance with an embodiment of the 
present invention; 

FIG. 3 is an example of an input snapshot showing radio frequency 
(RF) versus time of arrival, as received by a wideband receiver included in an 
embodiment of the present invention; 

FIG. 4 is an example of another input snapshot showing pulse width 
versus time of arrival, as received by a wideband receiver included in an embodiment 
of the present invention; 

FIG. 5 is yet another example of an input snapshot showing pulse 
amplitude versus time of arrival, as received by a wideband receiver included in an 
embodiment of the present invention; 
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FIG. 6 is still another input snapshot showing RF frequency versus 
pulse width, as received by a wideband receiver included in an embodiment of the 
present invention; 

FIG. 7 is a plot of RF frequency versus pulse width resulting after one 
iteration of clustering , as performed in accordance with an embodiment of the 
present invention; 

FIG. 8 is a plot of RF frequency versus pulse width resulting after four 
iterations of clustering, as performed by an embodiment of the present invention; 

FIG. 9 is a plot of RF frequency versus pulse width resulting after six 
iterations of clustering, as performed by an embodiment of the present invention; 
and 

FIG. 10 is a plot of RF frequency versus pulse width resulting after 
seven Iterations of clustering, as performed by an embodiment of the present 
Invention. 



DETAILED DESCRIPTION OF THE INVENTION 

The present invention provides an unsupervised iterative classification 
method for radar pattern recognition. The method is self-organizing and requires 
minimal input from human interaction. 

The method of the invention forms clusters from a set of input data 
(samples), where each cluster consists of very similar data (samples). The method 



ITDE-PAV103US - 6 - 

first defines a measure of pattern sinnilarity and establishes a rule for assigning 
individual samples to the domain of a specific cluster center. The invention uses the 
Euclidean distance between two data points x and z, 

D= ||x-z|| 

as a measure of pattern similarity. The smaller the distance, D, the greater is the 
similarity between x and z. 

After a measure of pattern similarity is selected, the method of the 
invention sorts (or partitions) samples into cluster domains. The Euclidean distance 
measure, D, lends itself to this procedure, because it is a good measure of proximity. 
However, because the proximity of two patterns is a relative measure of similarity, it 
is necessary for the invention to establish a threshold to define degrees of acceptable 
similarity for the clustering method. 

A performance-index is chosen to measure the degrees of similarity 
and a procedure is used which minimizes the chosen performance index. One such 
performance index is the sum of the squared errors resulting from the cluster, and is 
a proximity measure given by 

Nc 

7=1 xeSj 

Where Nc is the number of cluster domains (or simply clusters), Sj is the set of 
samples belonging to the jth domain, and 
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xeSJ 



is the sample mean vector of set Sj, with Nj representing the sample size of Sj. 

There are other performance indices used in the method of clustering 
the samples, such as: (1) average squared distances between samples in a cluster 
domain, (2) average squared distances between samples In different cluster 
domains, (3) indices based on a scatter matrix and (4) minimum and maximum 
variance indices. 



An embodiment of the invention will now be described based on an 
algorithm referred to as Iterative Self Organizing Data Analysis Techniques 
(ISODATA), which is disclosed by J. T. Tou and R. C. Gonzalez, Pattern Recognition 
Principles, Addison-Wesley, 1974, Chapter 3, pp. 75-109. The ISODATA algorithm, 
generally designated as 10, is also shown in FIG. 1, and is further described below. 

For a set of N samples, {xi, Xj x^ }, ISODATA clustering algorithm 

includes the following steps: 

Step 1: Specify various clustering parameters, as follows: 
K = number of cluster centers desired; 
On = the minimum number of samples allowed in a cluster; 



9s = standard deviation parameter; 
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Gc = lumping parameter; 

L = maximum number of pairs of cluster centers which may be lumped; 
I = number of iterations allowed. 

Step 2: Distribute the N samples among the present cluster centers, 
using the following relationship: 

X e Sj if II X - zj II < ||x - z| i = l,2,...,Nc; . i 9^ j 

for all X in the sample set. In this notation, Sj represents the subset of samples 
assigned to the cluster center zj. 

Step 3: Discard sample sets with fewer than Gn members. That is, if 
for any j, Nj < Gn, discard Sj and reduce Nc by 1. 

Step 4: Update each cluster center zj, j = l,2,...,Nc, by setting it equal 
to the sample mean of its members (Sj), as follows: 

where Nj is the number of samples in Sj. 

Step 5: Compute the average distance Dj of samples in cluster 
domain Sj from their corresponding cluster center, using the following relationship: 
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=^ZII>^-4 j = l,2....Nc 

■^^ j xe 

Step 6: Compute the overall average distance of the samples from 
their respective cluster centers, using the following relationship: 

1 Nc 

step 7: The following decisions are then made: 

(a) If this Is the last Iteration, set Gc = 0 and go to Step 11; 

(b) If Nc < K/2, then go to Step 8; 

(c) If this is an even-numbered iteration, or if Nc> 2K, go to Step 11 
otherwise continue. 

Step 8: Find the standard deviation vector oj = (aij, a2j, anj)' for 
each sample subset, using the following relationship: 




i = l,2,..n; j = l,2,...,Nc 



where n is the sample dimensionality, x.k is the ith component of the kth sample in 
Sj; Z|j is the ith component of Zj, and Nj is the number of sample in Sj. Each 
component of aj represents the standard deviation of the samples in Sj along a 
principal coordinate axis. 
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Step 9: Find the maximum component of each aj, j = l,2,...,Nc and 
denote it by aj^ax. 

Step 10: If for any oj^a^, j = 1,2, ...,Nc, there are oj^ax > Qs and 

(a) Dj > Dave and Nj > 2(0n + 1), 
or 

(b) Nc < K/2 

then split zj into two new cluster centers zj+ and zj", delete zj, and increase Nc by 1. 

Cluster center zj* is formed by adding a given quantity yj to the 
component zj which corresponds to the maximum component of csj, (aj^ax). Similarly, 
zj" is formed by subtracting from the same component of zj. One way of specifying 
yj is to let it be equal to a fraction of aj^ax, that is yj = kaj^ax with 0 < k < 1. 

If splitting took place in this step, then go to Step 2; otherwise 

continue. 

Step 11: Compute the pairwise distances Dij between all cluster 
centers, as follows: 



^u=lhi -^il i = l,2,...,Nc-l; j = i + l,...,Nc 
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Step 12: Compare the distances Dij against the parameter Gc. 
Arrange the L smallest distances which are less than Qc in ascending order, as 
follows: 



[DiJi.Di2j2, , DijL] 

where Diji <Di2j2 <....<DijL and L is the maximum number of pairs of cluster centers 
which may be lumped together. The lumping process is described below in Step 13. 

Step 13: With each distance Di|jk, there is associated a pair of cluster 
centers zlk and zjk. Starting with the smallest of these distances, perform a pairwise 
lumping operation, according to the following relationship: 

For k = 1, 2, L, if neither zik nor zjk has been used in lumping during 
this iteration, merge these two cluster centers, using the following 
relationship: 



Delete zi^ and zjk and reduce Nc by 1. 

It is noted that only pairwise lumping is allowed and that a lumped 
cluster center may be obtained by weighting each old cluster by the number of 
samples in its domain. It will be understood that since a cluster center can only be 
lumped once, this step may not always result in L lumped centers. 



ITDE-PAV103US 



- 12- 



Step 14: If this is the last iteration, the algorithm terminates. 
Otherwise, go to Step 1 If any of the process parameters requires changing at the 
user' s discretion, or go to Step 2 if the parameters are to remain the same for the 
next iteration. An iteration is counted every time the procedure returns to Step 1 or 



Based on a flowchart of the ISODATA algorithm, illustrated in FIG. 1, 
the inventors developed a computer program, as described below, for radar threat 
clustering and radar identification/recognition based on the clustering. The program 
was written in MATLAB, although other languages may have been used. 

A listing of the MATLAB program for clustering radar data samples is 
provided in the following tables. 

Table A, threat_gen_n.m, lists a program for generating a snapshot of 
the radars' pulse descriptive words (PDWs). The snapshot includes PDW mixes from 
multiple radar threats, as they may be intercepted by wideband receiver 21, as 
shown in RG. 2. Exemplary snapshots (80 ms duration) are shown in FIGS. 3-6 
(discussed later). 
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Table A. Program threat gen n.m 



dear all; 

colors=['r'/b'/k'/g'/m',V,'p']; 

% 

% 

% — Program threat_gen_n.m 
% 

% This program is used to generate threat signal PDW for isodata.m 

% 

% - 

% 

ss="; 
% 

% argument for threatlist of each row: 

% 

% freq,std_freq,aoa,std_aoa,pw,std_pw,pa,std_pa,toa,std_toa,pri,id 

% 

% 

threatlist=[900,200,35,10,8,2,-5,10,22,0,1120,l;... 

1050,100,40,10,14,2,-15,10,300,0,1820,2;... 
1150,150,45,10,4,2,-20,10,1700,0,1920,3;... 
1310,190,50,10,6.5,2,-35,10,700,0,1320,4;... 
1210,5,53,4,12,2,-45,10,1900,0,2920,5]; 

% 

%the weight set: wl=0.; w2=0.1; w3=0.1; w4=0.1; w5=0.1; w6=0.; 
% 



[row,col]=size(threatlist); 

for threatnum=l:row 

a=threat(threatlist(threatnum,:)); 

tmptxt=sprintfCthreat_%d=a/,threatnum); %dynamically name and 
variable 

eval(tmptxt); 

ss=cat(l,ss,a); 

end 

ss_sorted=sortrows(ss,l); 
sss= zeros(size(ss) ) ; 
sss=ss_sorted; 
[row,col]=size(sss); 
% 

% 

% 

% ploting input PDWs 
% 

% ss_sorted ( : ,1 ) - TOA 
% ss_sorted(:,2) - AOA 
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% ss_sorted ( : , 3 ) - RF Freq 
% ss_sorted(:,4) - Pulsewidth 
% ss_sorted(:,5) - Pulse Amp 
% ss_sorted(:,6) - Threat ID 

% 

% 

figure 

scatter(ss_sorted(:,4),ss_sorted(:,3),5,ss_sorted(:,6),'filled') 

titleC Input Threat PDWs in 80 ms Snapshot','FontSize',12,'FontWelqht','bold') 

xlabel('Pulsewidth (usee)'); 

ylabelCRF Frequency (MHz)'); 

% 

figure 

scatter(ss_sorted(:,l),ss_sorted(:,2),5,ss_sorted(:,6),'filled') 

titleC Input Threat PDWs in 80 ms Snapshot','FontSize',12,'FontWeighf,'bold') 

xlabel('Time of Arrival (usee)'); 

ylabel('Angle of Arrival (deg)'); 

% 

figure 

scatter(ss_sorted( : , l),ss_sorted( : ,3),5,ss_sorted( : ,6),'filled*) 

titleC Input Threat PDWs in 80 ms Snapshot','FontSize',12,'FontWeight'/bold') 

xlabeK'Timeof Arrival (usee)'); 

ylabel('RF Frequency (deg)'); 

% 

figure 

scatter(ss_sorted(:,l),ss_sorted(:,4),5,ss_sorted(:,6),'filied') 

titleC Input Threat PDWs in 80 ms Snapshot','FontSize',12,'FontWelght','bold') 

xIabelCTime of Arrival (usee)'); 

ylabelCPulsewidth (usee)'); 

% 

figure 

scatter(ss_sorted( : , l),ss_sorted( : ,5), 5,ss_sorted( : ,6),'filled') 

titleC Input Threat PDWs in 80 ms Snapshot','FontSlze',12,'FontWeight'/bold') 

XlabeK'Timeof Arrival (usee)'); 

ylabelCPuise Amplitude (dBm)'); 

% 

for colindex=2:(col-l) 

sigma=std(sss(:,colindex)); 

mean 1 =mean(sss( : ,eolindex)) ; 

sss(:,eolindex)=(sss(:,colindex)-meanl)./sigma*5; 
end 

% 

mean_value=[mean(ss(:,2)),mean(ss(:,3)),mean(ss(:,4)),mean(ss(:,5))]; 

std_value=[std(ss(:,2)),std(ss(:,3)),std(ss(:,4)),std(ss(:,5))l; 

n_total=length(sss); 

% 

dispC***** Input PDWs *****■) 
disp(ss_sorted) 

dispC***** Normalized PDWs **♦**') 
disp(sss) 

dispC***** Mean Value Vector ****♦•) 
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disp(mean_value) 

disp('***** STD Vector **♦**•) 

disp(std_value) 

disp('***** Sample Size *****•) 
disp(n_total) 

% 

% store generated PDW's 
% 

d lmwrite( 'test_pdw_n .out',sss) 
dlmwrite('mean_value_n.out',mean_value) 
dlmwrite('std_value_n.out',std_value) 
dlmwrite('n_total_n.out',n_total) 

% 
% 

% — end of program 



Each PDW, which is a vector, is composed of four components, 
describing an intercepted radar pulse, as follows: (1) time of intercept (or arrival), 
TOA, (2) radio frequency, RF, (3) pulse width, PW, and (4) pulse amplitude, PA. It 
will be appreciated that in other embodiments of the present invention, less or more 
than four components (dimensions) of each PDW may be selected. For example, 
other components may be pulse repetition interval (PRI), modulation type, etc. 



Referring to FIG. 2, system 20 includes wideband receiver 21 for 
receiving desired components of each radar 1-N. Also included Is processor 26, 
coupled to wideband receiver 21, for generating each PDW using PDW generator 22, 
normalizing each PDW using PDW normalizer 23 and clustering each normalized PDW 
into a respective cluster using ISODATA module 24. The ISODATA module executes 
the steps of the ISODATA algorithm. After a predetermined number of iterations of 
the ISODATA algorithm, the clusters of PDWs may be formed and provided to radar 
classifier module 25, which may be included in processor 26 or may be a separate 
module. By matching the clusters against stored table 27, the latter containing 
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identifications of known radars (threats and/or non-threats), the radars may be 
classified and identified. 

Each raw PDW is normalized by module 23 of system 20, using the 
following relationship: 

PDW^,, = [PDW^^ -PDW^^J / STDp^^ 

where PDWNor is the individual normalized PDW vector, PDW^aw is the individual PDW 
as intercepted by wideband receiver 21, PDWavg is the average PDW vector of the 
entire snapshot, and STDpow is the standard deviation vector calculated from PDWraw 
and PDWavc- 

Table B, threat.m, lists a MATLAB function called by threat_gen_n.m to 
generate the PDWs. 
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Table B. Program threat.m 



function a=threat(inputvector) 

% 

% 

% — Program threat.m 

% 

% This program is used to generate PDW for threat 
% 

% - 

% 

% Snapshot Size 

% 

nt=80000; 

% 

% parameters of each threat 
% 

freq = inputvector(l); 

std_freq=inputvector(2); 

aoa=inputvector(3); 

std_aoa=inputvector(4); 

pw=inputvector(5); 

std_pw=inputvector(6); 

pa=inputvector(7); 

std_pa = 1 nputvector(8) ; 

toa=inputvector(9); 

std_toa=inputvector(10); 

pri=inputvector(ll); 

ww=inputvector(12); 

% 

% generate threat PDW's with n types each with k samples 
% 

np2=floor((nt-toa(l))/pri); 

%disp(np2) 

% 

al=toa(l); 
for irun=l:np2 
al=al+pri; 

a2=aoa(l)+(rand(l)-0.5)*std_aoa(l); 

a3=freq(l)+(rand(l)-0.5)*std_freq(l); 

a4=pw(l)+(rand(l)-0.5)*std_pw(l); 

a5=pa(l)+(rand(l)-0.5)*std_pa(l); 

a6=ww(l); 

a(iain,:)=[al a2 a3 a4 a5 a6]; 
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Table 1, isodata_n.m, is the main program, wfiich executes Step 1, 
Step 7 and Step 14 of the ISODATA algorithm, executed by module 24 of FIG. 2. In 
the exemplary embodiment of the invention, the performance index for measuring 
similarity between two PDWs uses only two components, namely RF and PW. 

The Euclidean distance, Dij, between PDWs (PDWI and PDWj) is 
calculated as follows: 

Dy = w,(RFi -RFj)' + W2(pWi 

where (RF,, PW^and (RFj, PWj) represent PDWj and PDWj, respectively. Two weights, 
wi and W2 are used, as an example, to adjust the relative size of the cluster (or 
equivalently, the pairwise distance between cluster centers) to be generated In 
ISODATA. The relative size may be adjusted as a function of the overall frequency 
and pulse width deviations, which likely are related to the number of input radar 
threats of the input snapshot, or may be adjusted as a function of dedicated 
frequency bands in which advanced emitters may reside and need to be clustered 
into a specific cluster size. 

Referring to Table 1, six weights are listed (wi - Wg). All weights are 
set to zero, except W3 and W4, which are RF frequency and pulse width, respectively. 
It will also be appreciated that initially at the start of the ISODATA algorithm, the 
number of clusters may be assumed to be 1. Samples to far away from a center of 
this original cluster may then be dropped from the cluster and a new cluster may be 
formed from the dropped samples. 
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Table 1. Program ISODATA n.m 

% Program isodata_n.m 

% 

% Interactive Self-Organizing Data Analysis Technique Algorithm 
% (ISODATA) 

% 

% input Data to be clustered 

% 
% 

ss=dlmread('test_pdw_n.out'); % input pdws 

fmean=dlmread('mean_value_n.out'); % mean value vector 
fstd=dlmread('std_value_n.out'); % std value vector 
n_total=dlmread('n_totaLn.out'); % total number of input pdws 

% 
% 

np=n_total; 

ndim=6; 

% 

% parameters and initial conditions 

% 

wl=0.; %TOA 
w2=0.0; %AOA 
w3=0.1; %RF Freq 
w4=0.15; %Pulsewidth 
w5=0.; %PulseAmp 
w6=0.; %ID 
% 

% initial number of cluster = 1 

% 

nc=l; 
% 

z(l,:)=ss(l,:); 
% 

% 

% ***** Step 1 (specify parameters) ******* 

% - 

% 

k=18; % number of cluster centers desired 

theta_n=l; % min numbers allowed in a cluster 
% 

% - 

% 

% original settings: 

% 

%theta_s=0.75; % std parameter 
%theta_c=2; % lumping parameter 
% 
% 

theta_s=0.75; % std parameter 
theta_c=l; % lumping parameter 
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% 

% — 

% 

L=0; % max number of pairs of cluster centers allowed to be lumped 

Imax=12; % number of iterations allowed 

% 
% 

dispC***** SODATA Cluster Seeking Algorithm *****♦') 

% 

dispC***** Number of Data Points *♦****•) 
disp(np) 

% 

dispC***** Dimension of PDW **♦***•) 

disp(ndim) 

% 

dispC***** Input Normalized PDWs Data *****•) 
disp(ss) 

disp('***** Input Mean Value Vector *****•) 
disp(fmean) 

dispC***** Input STD Vector ****♦') 

disp(fstd) 

% 

dispC***** Number of Clusters Desired *»****•) 

disp(k) 

% 

dispC***** Min Number of Data Required in a Cluster *****♦') 

disp(theta_n) 

% 

dispC***** STD Parameter *****♦•) 

disp(theta_s) 

% 

dispC***** Lumping Parameter *****♦•) 

disp(theta_c) 

% 

dispC***** Max Number of Cluster-Pairs Allowed to be Lumped 

disp(L) ^ ' 

% 

dispC***** Max Number of Iterations Allowed ♦*****') 

disp(Imax) 

dispC •) 

dispC***** Start Isodata Algorithm ******•) 

disp('************************************'\ 

% 

% 

for irun =l:Imax o/^ number of total "do loops" to step 14 

% set initial condition for each loop 

% 

goto2=0; 
gotoll=0; 
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split=0; 
% 

% ^ 

% 

step2 
steps 
step4 
steps 
step6 
% 

% 

o/q ***** step 7 ************* 

% — 

% 

if (irun == Imax) 

theta_c=0; 

gotoll=l; 
end 
% 

if (irun Imax & (mod(irun,2) == o | nc >= 2*k)) 

gotoll=l; 
end 
% 

if (irun Imax & nc > k/2) 

gotoll=l; 
end 

% 

if gotoll ~= 1 

steps 

step9 

steplO 

nc=ncc; 

if split == 1 
goto2 =1; 

end 
end 
% 

if goto2 ~=1 

step 11 

step 12 
end 

end %main loop 

% 

% Summary of Clusters' Positions 

% 

disp(« 

dispC***** Summary of Clustering Results') 
disp('***********************************«) 
dispC •) 

dispC***** Clusters Positions ******) 
sc=5; 
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for ipr=l:nc 

z(ipr,2)=z(ipr,2)*fstd(l)/sc+fmean(l); 

z(ipr,3)=z(ipr,3)*fstd(2)/sc+fmean(2); 

z(ipr,4)=z(ipr,4)*fstd(3)/sc+fmean(3); 

z(ipr,5)=z(ipr,5)*fstd(4)/sc+fmean(4); 
end 
disp(z) 

disp('***** Clusters Members ****♦■) 
for ipr=l:np 

yy(ipr,2)=yy(ipr,2)*fstd(l)/sc+fmean(l); 

yy(ipr,3)=yy(ipr,3)*fstd(2)/sc+fmean(2); 

yy(ipr,4)=yy(ipr,4)*fstd(3)/sc+fmean(3); 

yy(ipr,5)=yy(ipr,5)*fstd(4)/sc+fmean(4); 
end 

disp(yy) 
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Table, 2, step2.m, lists the program to perform ISODATA Step 2. 



Table 2. Program STEP2.m 
% ***** Step 2 (cluster data based on Euclidean distances) ****** 

disp('***** Step 2 ****•) 
dispC*** Run Number ***') 
disp(irun) 
% 

% If there are more than 1 clusters 

% 

dispC*** Number of Clusters ****') 
disp(nc) 

> 1 % for multiple clusters, 

% the else at line 94 
% 

% Euclidean Distance 
% 

foridis=l:nc 
for kdis=l:np 

dis2=wl*(ss(kdis,l)-z(idis,l))^2 + w2*(ss(kdis,2)- z(idis,2))^2 + 
w3*(ss(kdis,3)-z(idis,3))^2 + w4*(ss(kdis,4)-z(idis,4))^2 + .. 
w5*(ss(kdis,5)-z(idis,5))^2 + w6*(ss(kdis,6)-z(idis,6))^2; 
dis(kdis,idis)=sqrt(dis2); 
end o/o for kdis at line 17 

end o/o for jdis at line 16 

% 

% Sorting into clisest cluster center 

% 

for kdis=l:np 

comp=dis(kdis,l); 
i_sort(kdis) = l; 
for idis=l:nc 

if dis(kdis,idis) < comp 
comp = dis(kdis,idis); 
i_sort(kdis)=idis; 
end o/o for jf at line 31 

end for jdis at line 30 

end o/o for kdis at line 27 

% 

% Count how many members in each cluster 
% 

for idis=l:nc 
n_total(idis)=0; 
for kdis=l:np 

if i_sort(kdis) == idis 

n_total(idis)=n_total(idis)+l; 
end o/o for if at line 43 
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% for kdis at line 42 
% for idis at line 40 



end 
end 

% 

% If there is only one cluster 
% 

else o/j, for if at line 12, only one cluster 

for kdis=l:np 
i_sort(kdis) = l; 

end o/o for kids at line 52 

n_total(l)=np; 
end o/o for else at line 51 

% 

% print present cluster centers 
% 

dispC******** Present Cluster Centers *»****■) 
% 

for iprint=l:nc 
z(iprint,:); 

end o/o for jphnt at line 62 

dlsp(z) 

% 

for iprint=l:np 

yy(iprint,:)=[ss(iprint,:) i_sort(iprint)]; 
end 

dispC ****** Present Cluster Members *******') 



dispC Input Data 

disp(yy) 
sc=5; 

for ipr=l:np 

yz(ipr,l)=yy(ipr,l); 

y2(ipr,2)=yy(ipr,2)*fstd(l)/sc+fmean(l) 

yz(ipr,3)=yy(ipr,3)*fstd(2)/sc+fmean(2) 

yz(ipr,4)=yy(ipr,4)*fstd(3)/sc+fmean(3) 

yz(ipr,5)=yy(ipr,5)*fstd(4)/sc+fmean(4) 

yz(ipr,6)=yy(ipr,6); 

yz(ipr,7)=yy(ipr,7); 

end 



Cluster') 



% plot clustering results after each run 
% 

% 

% 
% 
% 
% 
% 
% 
% 
% 
% 
% 



yz(:,l) - TOA 
yz(:,2) - AOA 
yz(:,3) - RF Freq 
yz(:,4) - Pulsewidth 
yz(:,5) - Pulse Amp 
yz(:,6) - Threat ID 
yz(:,7) - Cluster Number 
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% 

figure 

scatter(yz( : ,4),y2( : ,3),5,yz( : ,7)/filled') 
title(sprintf('Clustered PDWs, Run Number= 
5 %d',irun)/FontSi2e'42/FontWeight'/bold'); 
xlabel('Pulsewidth (usee)'); 
ylabel('RF Frequency (MHz)'); 
M(irun) = getframe; 



Table 3, step3.m, lists the program to perform ISODATA Step 3. 



10 Table 3. Program STEP3.m 

% ***** Step 3 (discard subsets with fewer than theta_n members) ****** 

dispC***** Step 3 ****•) 
n_new=0; 
15 foridis=l:nc 

if n_total(idis) >= theta_n 
n_new=n_new+l; 
n_total(n_new)=n_total(idis); 
end o/o for if at line 7 

20 end o/o for id is at line 6 

nc=n_new; 
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Table 4, step4.m, lists the program to perform ISODATA Step 4. 



Table 4. Program STEP4.m 
***** Step 4 (update cluster centers as their samples' means) ****** 

dispC***** Step 4 ****■) 
foridis=l:nc 

If n_total(idis) ~= 0 
z(idis,:)=0; 
forkdis=l:np 

if i_sort(kdis) == idis 

z(idis,:)=z(idis,:)+ss(kdis,:); 
end % for if at line 9 

end o/o for kdis at line 8 

z(ldis, : )=z(idis, : )/n_total(idis); 
end o/o for if at line 6 

end o/o for jdis at line 5 

% 

% print new cluster centers 

% 

dispC***** New Cluster Points **»***•) 
% 

%for iprint=l:nc 
% iprint, z(iprint,:) 

%end o/o for jphnt at line 21 

disp (z) 

% 
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Table 5, stepS.m, lists the program to perform ISODATA Step 5. 



Table 5. Program STEPS. m 

% ***** step 5 (compute average distance to cluster center within each 
% cluster) ****** 

% 

% 

dispC***** Step 5 ****•) 
for idis=l:nc 

if n_total(idis) > 0 
dis_c(idis)=0; 
for kdis=l:np 

if i_sort(kdis) == idis 

dls2=wl*(ss(kdis,l)-z(idis,l))^2 + w2*(ss(kdis,2)-z(idis,2))'^2 
w3*(ss(kdis,3)-z(idis,3))^2 + w4*(ss(kdis,4)-z(idis,4))^2 + . 
w5*(ss(kdis,5)-z(idis,5))'^2 + w6*(ss(kdis,6)-z(idis,6))'^2; 
dis_c(idis)=dis_c(idls)+sqrt(dis2); 
end % for if at line 10 

end o/o for kdis at line 9 

dis_c(idis)=dis_c(idis)/n_total(idis); 
end % for if at line 7 

end o/o for jdis at line 6 

% 

% print average distance within each cluster 
% 

dispC******* Ave Distance within each Cluster ******m 
% 

%foripr=l:nc 

% ipr, dis_c(lpr) 

%end 

disp(dis_c) 

% 
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Table 6, step6.m, lists the program to perform ISODATA Step 6. 



Table 6. Program STEP6.m 



o/o ***** sjgp g (compute overall average distance) ****** 
% 

dispC***** Step 6 ****') 

dis_t=0; 

foridis=l:nc 

if n_totai(idis) > 0 

dis_t=dis_t + n_total(idis)*dis_c(idis); 
end o/o for jf at line 8 

end % for idis at line 7 

% 

% print overall average distance 

% 

dispC******* Overall Ave Distance from Clusters ******'^ 

dis_ave=dis_t/np; 
disp(dis_ave) 

% 



(Table 7 Is skipped for convenience of numbering the tables so that 
they correspond to the numbered IDODATA steps.) 
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Table 8, steps. m, lists the program to perform ISODATA Step 8. 
Table 8. Program STEP8.m 

% 

% 

% ***** Step 8 (find STD vector for each sample set) ****** 

% ^ 

% 

dispC***** At Step 7, It decides to go to Step 8 *****•) 
% 

foridis=l:nc 
varl=0; 
var2=0; 
var3=0; 
var4=0; 
var5=0; 
var6=0; 
for kdis=l:np 

if i_sort(kdis) == id is 

varl=varl+wl*(ss(kdis,l)-z(idis,l))^2; 
var2=var2+w2*(ss(kdis,2)-z(idis,2))^2; 
var3=var3+w3*(ss(kdis,3)-z(idis,3))'^2; 
var4=var4+w4*(ss(kdis,4)-z(idis,4))'^2; 
var5=var5+w5*(ss(kdis,5)-z(idis,5))^2; 
var6=var6+w6*(ss(kdis,6)-z(idis,6))^2; 
end o/o for jf at line 16 

end % for if at line 15 

std(idis,l)= sqrt(varl/n_total(idis)); 
std(idis,2)= sqrt(var2/n_total(idis)); 
std(idis,3)= sqrt(var3/n_total(idis)); 
std(idis,4)= sqrt(var4/n_total(idis)); 
std(idis,5)= sqrt(var5/n_total(idis)); 
std (id is, 6) = sq rt(va r6/n_tota I ( idis) ) ; 
end o/o for idis at line 8 

% 

% print std vector for each cluster 
% 

' ****** STD Vector for each Cluster ******* 
% 

%for ipr=l:nc 
% ipr,std(ipr,:) 

"/"end o/o for jpr at line 37 

disp(std) 

% 
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Table 9, step9.m, lists the program to perform ISODATA Step 9. 

Table 9. Program STEP9.m 

% r : 

o/o ***** stgp 9 (fjpij ^f^g j^g^ component of each std vector and its 
orientation) ****** 

% - - 

% 

dispC***** step 9 *****') 
for idls=l:nc 

stdmax(idis)=std(idis,l); 

i_stdmax(idis) = l; 
end o/o foj. j^js at g 

for idis=l:nc 

if std(idis,2) > stdmax(idis) 

stdmax(idis)=std(idis,2) ; 

i_stdmax(ldis)=2; 
end 

if std(idis,3) > stdmax(idis) 

stdmax(idis)=std(idis,3); 

i_stdmax(idis)=3; 
end 

If std(idis,4) > stdmax(idis) 

stdmax(idis)=std(idis,4); 

i_stdmax(idis)=4; 
end 

if std(idis,5) > stdmax(idis) 

stdmax(idis)=std(idis,5); 

i_stdmax(idis)=5; 
end 

if std(idjs,6) > stdmax(idis) 

stdmax(idis)=std(idis,6); 

i_stdmax(idis)=6; 
end 

end o/o for jdis at line 10 

% 

% print max std vector and direction of each cluster 

% 

dispc****** Max STD and Its Direction *******) 
for ipr=l:nc 

max_std(ipr,:) = [lpr,stdmax(ipr),i_stdmax(ipr)]; 
end o/o for jpr at line 35 

disp(max_std) 
% 



ITDE-PAV103US 



- 31 - 



Table 10, steplO.m, lists the program to perform ISODATA Step 10. 



Table 10. Program STEPlO.m 



% 



o/o ***** step 10 (determine whether to split into more clusters) ****** 

5 % 

% 

dispC***** Step 10 *****•) 
% 

dispC******* New Splitted Clusters ******•) 
10 ncc=nc; 

foridis=l:nc 

if stdmax(idis) > theta_s 

if ((dis_c(idis) > dis_ave & n_total(idis) > 2*(theta_n+l)) | (nc < k/2)) 
gamma=0.5*stdmax(idis); 
15 split=l; 

if i_stdmax(idis) ==1 

z(idis,l)=z(idis,l)+gamma; 
z(ncc+l,l)=z(idis,l)-2*gamma; 
z(ncc+l,2)=z(idis,2); 
20 z(ncc+l,3)=z(idis,3); 

z(ncc+l,4)=z(idis,4); 
z(ncc+l,5)=z(idis,5); 
z(ncc+l,6)=z(idis,6); 
disp(z(ldis,:)) 
25 disp(z(ncc+l,:)) 

ncc=ncc+l; 

elseif i_stdmax(idls) == 2 

z(idis,2)=z(idis,2)+gamma; 
30 z(ncc+l,2)=z(idis,2)-2*gamma; 

z(ncc+l,l)=z(idis,l); 

z(ncc+l,3)=z(idis,3); 

2(ncc+l,4)=z{idis,4); 

z(ncc+l,5)=z(idis,5); 
35 z(ncc+l,6)=z(idis,6); 

disp(z(idis,:)) 

disp(z(ncc+l,:)) 

ncc=ncc+l; 

40 elseif i_stdmax(idis) == 3 

z(idis,3)=z(idis,3)+gamma; 

2(ncc+l,3)=z(idis,3)-2*gamma; 

z(ncc+l,l)=z(idis,l); 

z(ncc+l,2)=z(idis,2); 
45 2(ncc+l,4)=z(idis,4); 

z{ncc+l,5)=z(idis,5); 

z(ncc+l,6)=z(idis,6); 
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disp(z(idis,:)) 

disp(z(ncc+l,:)) 

ncc=ncc+l; 

elseif Lstdmax(idis) ==4 
2(idis,4)=z(idis,4)+gamma; 
z(ncc+l,4)=z(idis,4)-2*gamma; 
z(ncc+l,l)=z(idis,l); 
z(ncc+l,2)=z{idis,2); 
z(ncc+l,3)=z(idis,3); 
z(ncc+l,5)=z(idis,5); 
z(ncc+l,6)=z(idis,6); 
disp(z(idis,:)) 
disp(z(ncc+l,:)) 
ncc=ncc+l; 

elseif l_stdmax(idis) ==5 
z(jdis,5)=z(idis,5)+gamma; 
z(ncc+l,5)=z(idis,5)-2*gamma; 
z(ncc+l,l)=z(ldis,l); 
z(ncc+l,2)=z(idis,2); 
z(ncc+l,3)=z(idis,3); 
z(ncc+l,4)=z(ldis,4); 
z(ncc+l,6)=z(idis,6); 
disp(z(idis,:)) 
disp(z(ncc+l,:)) 
ncc=ncc+l; 

else 

z(idis,6)=z(idis,6)+gamma; 

z(ncc+l,6)=z(idis,6)-2*gamma; 

z(ncc+l,l)=z(idis,l); 

z(ncc+l,2)=z(idis,2); 

z(ncc+l,3)=z(idis,3); 

z(ncc+l,4)=z(idis,4); 

z(ncc+l,5)=z(idis,5); 

disp(z(idis,:)) 

disp(z(ncc+l,:)) 

ncc=ncc+l; 

% for if at line 14 
% for if at line 12 
% for if at line 11 
% for idis at line 10 



end 
end 
end 
end 
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Table 11, stepll.m, lists the program to perform ISODATA Step 11. 
Table 11. Program STEPll.m 



o/o ***** step 11 (compute pairwise distance between cluster centers and whether 
% ***** to lump) ****** 

% , 

% 

dispC***** Step 11 *****') 
for idisl = l:nc-l 
for idis2=2:nc 
if idis2 > idisl 

diss=wl*(z(idisl,l)-z(idis2,l))^2 + w2*(z(ldisl,2)-z(idis2,2))^2 + . 
w3*(z(idisl,3)-z(idis2,3))^2 + w4*(z(idisl,4)-z(idis2,4))^2 + 
w5*{z(idisl,5)-z(idis2,5))^2 + w6*(z(idisl,6)-z(idis2,6))A2; 
dis_bw(idisl,idis2)=sqrt(diss); 
Iump(idisl,idis2)=0; 
if dis_bw(idisl,idis2) <= theta_c 

Iump(idisl,idis2) = l; 
end o/q for if at line 15 

end % for if at line 9 

end o/o for jdjs2 at line 8 

end o/o for jdisl at line 7 



ITDE-PAV103US 



34 



Table 12, stepl2.m, lists the program to perform ISODATA Step 12 



and Step 13. 



Table 12. Program STEP12.m 



% 



o/o ***** step 12 & 13 (lump L or less clusters as allowed and updating new 
o/o ***** cluster centers) ****** 

% - - 

% 

dispC***** Step 12 and 13 *****•) 
for idisl=l:nc-l 
for idis2=2:nc 

if (idisl ~= idls2 & idis2 >= idisl & Iump(idisl,idis2) == 1) 

% 

zl_n=n_total(idisl)*z(idisl,l)+n_total(idis2)*2(idls2,l) 
z(idisl,l)=zl_n/(n_total(idisl)+n_total(idis2)); 
% 

z2_n=n_total(idisl)*z(idisl,2)+n_total(idis2)*z(idis2,2) 
z(idisl,2)=z2_n/(n_total(idisl)+n_total(idis2)); 

% 

z3_n=n_total(idisl)*z(idisl,3)+n_total(idis2)*z(idis2,3) 
z(idisl,3)=z3_n/(n_total(idisl)+n_total(idis2)); 
% 

z4_n=n_total(idisl)*z(idisl,4)+n_total(idis2)*z(idis2,4) 
z(idisl,4)=z4_n/(n_total(idisl)+n_total(idis2)); 
% 

z5_n=n_total(idlsl)*z(idisl,5)+n_total(idis2)*z(idls2,5) 
z(idisl,5)=z5_n/(n_total(idisl)+n_total(idis2)); 

% 

z6_n=n_total(idisl)*z(idisl,6)+n_total(idis2)*z(idis2,6) 
z(idisl,6)=z6_n/(n_total(idisl)+n_total(idis2)); 
% 

disp(z(ldisl,:)) 
2{idis2,:)=99999; 
end o/o for jf at line 9 

end o/o for icjis2 at line 8 

end o/o for idisl at line 7 



To illustrate the operation of the invention, a simple test case of an 
electronic warfare (EW) scenario consisting of five (5) radar threats was provided to 
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a simulation of system 20. Tlie five radar threats and their characteristics are listed 
in Table 13. 



Table 13. Radar Threat Characteristics for ISODATA 
Threat RF Freq PRT 


Simulation 
Pulse Width 


Pulse Amplitude 


1 


800 - 1,000 MHz 


1120 us 


7 - 9 us 


0 to -10 dBm 


2 


1000 - 1100 MHz 


1820 us 


13 - 15 us 


-10 to -20 dBm 


3 


1075 - 1225 MHz 


1920 us 


3 -5 us 


-15 to -25 dBm 


4 


1215 - 1405 MHz 


1320 us 


5.5 - 7.5 us 


-30 to -40 dBm 


5 


1207.5 - 1212.5 MHz 


2920 us 


11 -13 us 


-40 to -50 dBm 



Snapshots of PDWs of this EW scenario were generated by the program 
threat generator listed in Table A. FIG. 3 illustrates a snapshot of the radars' RF 
frequency (RF) versus TOA. FIG. 4 illustrates a snapshot of the radars' puise width 
(PW) versus TOA. FIG. 5 illustrates a snapshot of the radars' pulse amplitude (PA) 
versus TOA. These figures indicate that radar pulses from multiple threats are 
overlapped in RF, PW and PA. There, indeed, are five (5) radar threats in the input 
snapshots, when the snapshots are shown plotted as RF versus PW, as illustrated in 
FIG. 6. 

Performance of system 20 in clustering and classifying the five radar 
threats is summarized in HGS. 7-10. These figures illustrate the result and progress 
of clustering after 1, 4, 6 and 7 iterations, respectively. The system performs well 
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and after seven iterations, the system clusters the input snapshots into five radar 
threats, without errors, as shown in FIG. 8. 

It will be appreciated that system 20 may be used to cluster EW 
scenarios consisting of mixes of stable radars and advanced radars, such as dwell 
switched and frequency agile radars. 

It will also be appreciated that to cluster advanced emitters having 
frequency agility capability, the weighs (wi, wj, and others, if necessary) used in 
Euclidean distance calculations between PDWs may be made adaptive, so that PDWs 
from different threats may be sorted into different clusters and PDWs from the same 
threat will not be partitioned into multiple clusters. As an example, the weights may 
be made a function of the operational frequency band of the radar emitter and the 
size of clusters generated may be adjusted to prevent threat splitting. 

Although illustrated and described herein with reference to certain 
specific embodiments, the present invention is nevertheless not intended to be 
limited to the details shown. Rather, various modifications may be made in the 
details within the scope and range of equivalents of the claims and without departing 
from the spirit of the invention. 



